Optimal basis set for electronic structure calculations in periodic systems 
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An efficient method for calculating the electronic structure of systems that need a very fine sampling 
of the Brillouin zone is presented. The method is based on the variational optimization of a single 
(i.e. common to all points in the Brillouin zone) basis set for the expansion of the electronic orbitals. 
Considerations from k -p-approximation theory help to understand the efficiency of the method. The 
accuracy and the convergence properties of the method as a function of the optimal basis set size 
are analysed for a test calculation on a 16-atom Na supercell. 



I. INTRODUCTION 

First-principles methods have become a widespread 
tool for studying the microscopic aspects of a wide class of 
condensed systems, including their crystalline phases. No 
empirical modelling of the interatomic interactions is re- 
quired in first-principles simulations, since the electronic 
ground state is evaluated at each atomic configuration!!!. 
In a large set of systems, mostly in metals, but also 
in some molecular crystals like hydrogeno, a very fine 
sampling of the Brillouin Zone (BZ) is required to pro- 
vide a sufficiently accurate description of the electronic 
ground state. Since the computational load for refin- 
ing BZ sums grows linearly with the number of sam- 
pling points (k-points), special k-point grids that min- 
imize sampling errors, as well as interpolation techniques 
(the so called tetrahedron methodcl) have been intro- 
duced in the pastoQ. However, little attention has been 
paid to the fact that Bloch states at nearby k-points are 
not totally uncorrelated. In fact, once the Bloch states 
at a given k-point are known, those at a closeby k-point 
can be easily estimated by perturbation theory, within 
the so-called k ■ p approximation!! This extremely help- 
ful piece of information has been little exploited so far 
in standard first- principles approaches, where for each 
k-point a full solution of the Kohn-Sham equations is 
instead calculated independently, e.g. through diagonal- 
ization of the Kohn-Sham hamiltonian on a finite basis 
set (typically plane waves in the case of crystalline sys- 
tems) . An approximate technique based on the exact so- 
lution of the KohrtSham hamiltonian at a few k-points 
has been proposecH, but the analysis of the errors asso- 
ciated with such method, as well as their correction, has 
proved rather difficult!]. 

In this work we propose a simple and efficient method 
for solving the Kohn-Sham equations for a very large 
number of k-points in the BZ, with an effort that in most 
cases is comparable to that of solving the problem for 
just a few k-points. The method is based on a rigorous 
formulation of the ideas sketched above, and relies on 
the construction of a single, optimal basis set for the 
expansion of the Bloch orbitals, which is common to all 
k-points. The application of this method to the study of 
the compressed phases of molecular hydrogen has been 



recently reportedEErO. In Section II we formulate the 
method. In Section III we analyse the reliability of the 
method on a test case (16 Na atoms in the bcc lattice), 
and in Section IV we analyse the computational cost. 
Section V contains the conclusions. 



II. THE METHOD 

In first-principles simulations the electronic ground 
state is calculated within the Kohn-Sham, self-consisteet 
one-electron formalism of density functional theoryll3. 
Translational invariance is introduced by averaging the 
electronic density and the total energy of the system on 
a finite number of k-points in the BZ: 
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P (r)=5> k £ /i(k) | yjf(r) | 2 , (1) 
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E[p]=T + E cxt \p] + E H [p] + Exc[p] (2) 

with 

T = f>k /*( k ) / # k ( r ) ^ k ( r ) dr - 

k i " ^ ' 

(3) 

where i labels bands, k a k-vector in the BZ, lu^ its 
weight, /i(k) are occupation numbers, and ipf (r) are the 
Kohn-Sham orbitals. The last three terms in (|J) are 
the external, Hartree, and exchange-correlation contri- 
butions to the energy, and depend on the BZ averaged 
electronic density (|lj). Atomic units are used throughout 
this work. 

Variational equations for the periodic part of the 
Kohn-Sham orbitals uf (r) = exp (— ik • r)tpf(r) are read- 
ily obtained from expression (|2|): 

(-¥ - zk • V + ^ + V\p(r)]j uftr) = efuf(r) , 

(4) 
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where V[p] is the Kohn-Sham potential, containing ex- 
ternal, Hartree and exchange-correlation contributions. 
Normally, equations (||) are solved independently for 
each k-point, a new density is constructed via (Q), and 
the cycle is iterated to self-consistency. Equations (^) 
are often solved by expanding the orbitals correspond- 
ing to the various k-points in a basis set of plane waves, 
which remains unchanged when the atomic configuration 
is changed. The choice of a plane- wave basis set offers 
some practical advantages, including the convenience of 
calculating a number of quantities (matrix elements, form 
factors, etc) a single time throughout the calculation. 
However, working with a plane-wave basis set becomes 
particularly demanding from the computational point of 
view when a large supercell and/or a large number of 
fc-points has to be considered, a situation which is of- 
ten encountered when dealing with metals. Moreover, a 
plane-wave basis set cannot straightforwardly exploit the 
similarity of Bloch states at closeby k-points. 

In the following we will show that relaxing the require- 
ment of a "simple" but inmutable basis set in favor of a 
variational search for the "best" basis set at each atomic 
configuration may improve substantially the efficiency of 
BZ summations. To this aim, let us expand the Kohn- 
Sham orbitals on a generic basis set fa (i = 1, ...,N) of 
size N: 
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^(r) = V Mk)^-(r) 



(5) 



In the case of standard calculations fa would be a plane 
wave - 4>i(r) ~ exp(ig • r), being g a reciprocal lattice 
vector - and solving the Kohn-Sham equations would 
amount to minimize the energy functional (^|) with re- 
spect to the expansion coefficients ay, without modify- 
ing the basis set. In the present approach we minimize 
the energy functional (^|) with respect to both the expan- 
sion coefficients {(%■} and the basis set functions {4>i}. In 
other words, we look simultaneously for the ground state 
orbitals and for basis set that better describes them. For 
the sake of clarity, we note that this would be a redun- 
dant request for calculations with a single k-point, where 
a single set of Kohn-Sham equations have to be solved. 
The request of an optimal basis set becomes non redun- 
dant when the basis set is used to describe Kohn-Sham 
equations at many k-points. 

Explicit minimization of (|^) with respect to, simulta- 
neously, a,ij and fa leads to the following system of equa- 
tions. For the expansion coefficients we have the usual 
eigenvalue equations: 
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For the basis-set functions we have 
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N 



&,,(k) = E/*(k) 4(k)a a (k) 



(11) 



The A coefficients in @ are Lagrange multipliers in- 
troducedj-to enforce the orthonomalization of the basis 
functionstl. The electron density in ^ is constructed 
from the Kohn-Sham orbitals (||) using ([y]), as 
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where 



&jl = E Wk M k ) 



(12) 
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More detailed expressions for the above equations as well 
as for atomic forces and stress are given in the Appendix, 
for the case of pseudopotential plane wave (PPW) calcu- 
lations. 

The problem of simultaneously solving Kohn-Sham 
equations for k-points independently, as done in stan- 
dard electronic structure calculations, is here replaced by 
the problem of solving a single set of Kohn-Sham-like 
equations, eqs. (|^), supplemented with Nk diagonaliza- 
tions of the N x N matrix in (^) . As we will show in the 
next Section, the number of basis functions N required 
for the optimal basis set {fa} can be chosen to be much 
smaller than that in the basis set used to describe each 
single fa, e. g. the number of plane wave components in 
PPW calculations. This implies that, if Nk is large, the 
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computational effort required by this procedure is much 
smaller than the effort required to solve the standard set 
of Kohn-Sham equations at all k-points. The main lim- 
itation of the method resides, then, in the size N of the 
optimal basis set. If N turns out to be comparable to 
the original number of basis functions, then there is no 
point on pursuing this strategy. 

While in standard calculations the number of Kohn- 
Sham equations is roughly (exactly for insulators) half 
the number of electrons, in our approach N has to be suf- 
ficiently large to guarantee the proper description of the 
Kohn-Sham orbitals (|5|) at all k-points. A full analysis 
of the convergence of the method as a function of N will 
be presented later, for a specific case. There are however 
reasons to expect that N need not be too large. In fact, a 
reasonable approximation to the exact optimal basis set is 
given by the Kohn-Sham orbitals at one selected k-point, 
say at zone center (T). We call this basis set: {<j>f }. Since 
the optimal basis set is defined as the best basis set in a 
variational fashion, {0^} will necessarily provide a worse 
description of the Kohn-Sham orbitals (except, of course, 
at r) than {<j>i}. The convergence properties of {</>j} will 
thus be better than those of {<p\ }. If the BZ of the sys- 
tem of interest is sufficiently small (which is equivalent 
to ask that the unit cell of the calculation is large), then 
k ■ p perturbation theory suggests that Kohn-Sham or- 
bitals at any k-point other than T may be expressed in 
terms of a few number of energetically close T-point or- 
bitals {<j)f }Q. In other words, only a few excited T-point 
orbitals are required for a reasonable description of the 
orbitals at all other k-points. This property has been 
already noted and used in the ocmiepct of first-principles 
electronic structure calculationsfflli-J. Xiur approach dif- 
fers from that of Robertson and Paynetm in the fact that 
their basis set is given by the Kohn-Sham orbitals at T 
(or at a small set of k-points), while our method includes 
a variational search for the best basis set. A comparison 
of the convergence properties of the two approaches is 
given in the next Section. 

We also observe that our method has only one vari- 
ational parameter, namely the size N of the basis set. 
This implies that systematic improvements of the accu- 
racy can be achieved by simply increasing N. 

Finally, we notice that eqs. (^), like ordinary Kohn- 
Sham equations, have to be solved self-consistently, since 
the Kohn-Sham potential Vks depends on the density, 
and thus on wavefunctions. However, unlike ordinary 
Kohn-Sham equations, eqs. ([)]) depend also on the ex- 
pansion coefficients ay. The ordinary self-consistent it- 
erative procedure has, consequently, been modified in the 
following way: 

1. an initial guess of the density is provided, e.g. as 
a sum of atomic densities; the ay coefficients are 
initially set to <5y , 

2. eqs. ([}]) are solved, 



3. the orbitals (which are no longer Kohn-Sham or- 
bitals) are used to construct A°- and py, 

4. eq. @ is solved for each k-point, via direct diago- 
nalization, 

5. a new density is constructed, 

6. back to 2. 

Notice that this procedure does not involve a major 
implementation effort, but only a few modifications to 
the usual self-consistency methods. 

III. CONVERGENCE OF THE OPTIMAL BASIS 

SET 

In order to analyse the convergence properties of the 
method we have chosen the case of Na metal, whichjs 
properly described with a soft, local pseudopotentialEJ. 
The Ceperley- Alder local density functional was adopted 
for the exchange-correlation termlill Sixteen Na atoms 
in the bee structure have been placed in a simple cubic 
supercell (2x2x2 bec conventional cells), and the full 
Brilloum. zone of the supercell has been sampled with 108 
k-pointso. A plane wave energy cutoff of 15 Ry was cho- 
sen, amounting to 1596 plane waves in the expansion of 
the optimal basis set. The convergence of the method has 
been analysed by performing calculations with N = 19, 
32, 51, 81, and 1596EZJ. Since our aim is at checking the 
convergence properties of the method, and not the over- 
all accuracy as determined by the other approximations 
(LDA, k-point grid, etc.), we will consider as converged 
the values obtained with N = 1596 and Nk = 108. As 
a comparison, we have also determined the convergence 
properties of a calculation where the basis set, instead of 
being "optimal" , is chosen to be the set of Kohn-Sham 
orbitals resulting from a self-consistent calculation using 
the |Tppoint only. This correponds to the approach of 
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FIG. 1. Convergence of the total energy of solid bcc Na 
using a 16-atom supercell, as a function of the number of the 
basifiiiset size N. Dashed lines correspond to the approach of 
Reffl'H, and solid lines to the present method. The horizontal 
line is the converged value. 

In Figs. 1 to 3 we show the convergence of the different 
physical properties as a function of N, namely the total 
energy, stress, and electronic density of states. The total 
energy (Fig. 1) shows an extremely good convergence 
even for the smallest value of N. The accuracy is better 
than 0.5 % for N = 19, and better than 0.1 % for N = 32 
states. The same behavior is observed in the force con- 
stant calculated for a distortion along an optical phonon: 
the force is already converged to better than 0.1 % using 
32 states in the expansion. Thus, the errors introduced 
in energies and forces due to truncation of the expansion 
are much smaller than, e.g., the typical accuracy of the 
usual density functionals used, namely the LDA or the 
GGA. Stresses also show good convergence properties, as 
reported in Fig. 2. Variations of stress (Act) with respect 
to N should in this case be compared with the value of 
the bulk modulus of Na (B = 64 kbar). In fact, an error 
Ac in the stress yields an error in the equilibrium lat- 
tice spacing that can be estimated as Aa/a ~ Acr/3-B. 
Thus, even with N = 19 the equilibrium lattice spacing 
is converged within less than 1 % . 
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FIG. 2. Convergence of the trace of the stress tensor of 
solid bcc Na using a 16-atom supercell, as a function of the 
basiBi|Set size N. Dashed lines correspond to the approach of 
Ref.u'u, and solid lines to the present method. The horizontal 
line is the converged value. 

The electronic density of states (EDOS) is also suf- 
ficiently accurate for the relevant portion of the spec- 
trum, namely the one below the Fermi energy (Ep), us- 
ing iV = 32. In contrast, N — 19 appears to give an 
inaccurate EDOS, particularly in the vicinity of Ep, as 
expected in a k-p picture. In fact, as discussed in Section 



II, a small basis set is more likely to affect the accuracy in 
the description of states with higher energy, i.e. closer to 
Ep, due to the insufficient number of high energy states 
in the optimal basis set. 
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FIG. 3. Convergence of the electronic density of states of 
bcc Na using a 16-atom supercell, as a function of the ba- 
sis set size N. The different curves correspond to different 
values of N: long-dashed (N = 19), short-dashed (N = 32), 
dot-dashed (N = 51), and solid (TV = 81). Only results ob- 
tained with the present method are reported. 

We now compare the results obtained with our_method 
against those obtained with the method of Ref.Era. Our 
energies are always lower, for a given basis set size JY., 
than energies calculated using the method of Ref.Em. 
This is not unexpected, since the optimal basis set is 
variationally the besl basis set. Moreover, we find that 
the method of RcfBEl provides results with an accuracy 
similar to ours only when the size of the basis set is 30-50 
% larger than ours. The most striking difference, how- 
ever, is in the convergence of the stress tensor. While our 
method converges very smoothly to the exact value, the 
method of Ref.ErEl shows an erratic behavior, that clearly 
prevents every attempt to determine the exact value on 
the basis of size scaling. We believe that such a marked 
difference in the convergence properties has to be traced 
to the variational nature of our approach. Therefore, our 
method not only improves on that of Ref.Blj in the rate 
of convergence (which, in the case of energy, results in a 
30-50 % saving in computational time, as shown in the 
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next Section) , but shows a much smoother behavior as a 
function of N, resulting in a more accurate estimate of 
the error introduced by the finite size of the basis set. 

IV. ESTIMATE OF THE COMPUTATIONAL 
LOAD 

The computational effort in a standard plane wave 
calculation scales as NkNf, an( isNp W , where Nb an d s is the 
number of occupied bands (equal or slightly larger than 
half the number electrons) and N pw is the size of the 
plane wave basis set. This estimate is based on the 
assumption that standard diagonalization methods are 
used, such as iterative diagonalization, where the effort of 
determining the first N^ands eigenvalues of a N pw x Npw 
matrix scales like Nb an d s Np W . Since every k-point re- 
quires a separate diagonalization, the total effort scales 
as N k N bands N^ w . 

In the present approach, if we assume that the effort 
of diagonalizing Nk times the N x N matrix (||) is negli- 
gible with respect to the effort required to solve eqs. 
then the load scales as NNp W . Thus, a gain of a fac- 
tor NkNbands/N is expected with our method. Since in 
typical applications N ~ 2 4 Nb an ds, and Nk ~ 10 2 , a 
computational gain of more than one order of magnitude 
is obtained. The neglected term becomes increasingly im- 
portant as the size of the system (N) increases, scaling as 
Nk x iV 3 . Therefore, it becomes dominant for sufficiently 
large N, thus discouraging the use of the present method- 
ology in the case of large unit cells and BZs sampled in 
a few k-points. 

The above estimate is based on the assumption that 
solving eqs. @ is as time consuming as solving the stan- 
dard Kohn-Sham equations, except for the fact that the 
requested number of states is larger. Ohs experience on 
Na (this work) and molecular hydrogerJiSii 2 ] shows that 
the computational overload of solving eqs. (g) basically 
scales linearly with the number of states; the extra ef- 
fort required for solving a single state is negligible when 
compared to that of solving one of the usual Kohn-Sham 
states. Nevertheless, we are currently not in a position 
to generalize this conclusion to a generic system. 

V. SUMMARY AND CONCLUSIONS 

We have presented and analysed the convergence prop- 
erties of a method for performing electronic structure cal- 
culations and first-principles molecular dynamics simula- 
tions for systems that need a very fine sampling of the 
electronic Brillouin zone. The method is based on the 
construction of a basis set which optimally describes the 
Kohn-Sham orbitals at all k-points. The rapid conver- 
gence of the method on the size N of the basis set is 
connected with the properties of the k ■ p approximation 
and has been fully analysed in the case of a sixteen Na 



atoms supercell. We have shown that N can be kept rea- 
sonably small, of the order of a few times the number 
of occupied bands (2 to 4 times, depending on the accu- 
racy desired), resulting in a gain in computational load of 
about one order of magnitude, with respect to standard 
methods. 

Even if our analysis was mainly done in connection 
with PPW calculations, the same philosophy can be 
adopted in othec electronic structure schemes like the 
LAPW methodlia, provided that the number of orbitals 
required in the optimal basis set remains significantly 
smaller than the number of basis functions in the origi- 
nal basis set. 

In conclusion, this technique constitutes a useful tool 
for studying the electronic, structural and dynamical 
properties of metals and more generally for any system 
whose proper description requires a fine sampling of the 
BZ, like, for example, molecular hydrogen. 
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VII. APPENDIX 

In this appendix we give some useful expressions for 
the total energy, atomic forces and stress, in the case of 
pseudopotential plane wave calculations. 

The different contributions to the total energy (g) can 
be divided into the class of terms that only depend on 
the density (the hartree, the exchange-correlation, and 
the local part of the pseudopotential), and the class of 
terms that explicitly depend on the wavefunctions (the 
kinetic and the nonlocal pseudopotential contributions). 
While the former can be calculated in the usual way once 
the density is known (see below for an efficient way to 
evaluate density), the latter can be expressed more con- 
veniently as follows. For the kinetic energy we have 

N N f 2 . 

r = E"*EE 6 ti 4 M + k ■ pv + Y 5ij = 

k i=l j=l J 
N 

= E {b^M + b^} + 

8=1 

N N 

+ EE {bS^M-< } -^} - ( 14 ) 

with 

I ~|&> , (15) 
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and for the non-local contribution we have 
a i 

4l=EE a °™ E E^Mk) 

s—l m— — l ij k 



(16) 



1=1 



where, in a plane waves basis set, _F/jz m (k) assumes the 
following form: 



FljlmW — 



,(g + k) exp(ig-R 7 ) &(g) , (17) 



being {g} the reciprocal lattice vectors and Rj the coor- 
dinates of atom / belonging to species s. The coefficient 
is the projection of the g-th plane wave component 
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expanded around the k- vector k, onto-thc ^m-th subspace 
according to the Klcinman-BylandeitJ prescription: 
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(g + k) =< 8vt<p s lm Y lm | g + k > 



(18) 



with Svf the Z-th angular momentum pseudopotential 
component for species s, tp s lm the corresponding pseu- 
doatomic orbital, and Yi m the spherical harmonic func- 
tions. Equation ( |ilf ) also contains a normalization fac- 
tor, which is given by the expression af m —< 5vftpf m | 
(pf > _1 . N s is the number of atoms of species s, and 
a is the number of species described through angular- 
dependent pseudopotentials. 

It is interesting to remark that, since now the factors 
w im depend explicitly on the k-vector, it is no longer pos- 
sible to perform the BZ summation on the coefficients 
fr«(k) once and for all, as it is done for the kinetic and 
local potential terms. Now this sum has to be carried 
out for each PW component g, and this is likely to re- 
sult in a costly computational scheme. So far we have 
not studied this issue in detail, because this method has 
only been applied to hydrogen and sodium, where a local 
pseudopotential description is feasible, but we hope that 
alternatives can be found so that the non-local part does 
not become the bottleneck of the calculation. 

The calculation of the density is a subtle issue, because 
expression (|l^) involves a double sum of matrix products 
over all states, occupied and empty. Even if it has been 
suggestcdLJ that this operation is computationally conve- 
nient, we have found it more efficient to first transform 
the orbitals at T by multiplying them with the rotation 
matrices a,j(k), and then carry out a single summation 
over the occupied states. 

The last term in (|) is easier to compute in reciprocal 
space, where it assumes the following expression: 

(g + k)J^ exp(ig-R / ) F Ijlm (k) . 
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The forces on the nuclei are trivially unchanged from 
the original ones, as long as local pseudopotentials are 
used. This is because the electronic kinetic term does 
not depend explicitly on the nuclear positions. In the 
case of angular-dependent pseudopotentials the contri- 
bution from the last term in ( |l6|) to the force on atom 
/ belonging to species s, is modified from the standard 
result in the following way: 

2 E a "™ E E Wk x 



dRr 



(20) 



The calculation of the stress matrix is also modified 
with respect to standard calculations only for the kinetic 
(kin) and nonlocal pseudopotential (NL) contributions. 
The resulting expression for the kinetic contribution is: 



'a/3 



with 



and 
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^£ = E^^^(k) (23) 

k 

B ia/3 = E Wk ka k P b "(k) 
k 

The nonlocal contribution is instead: 

a 

- 2 EE^» E E wk M k ) x 

ij k 



NL 
T aj3 



s=l lr, 



N. 



ERe( ^Jk)^ m(k; 



1=1 



dh af: 



(24) 



where is the matrix of the three primitive Bravais 
vectpjrsEII, and the partial derivative in (E4|) is given in 
Ref.El 



i=i 
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